Correlation energy and spin polarization in the 2D electron gas 
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The ground state energy of the two-dimensional uniform electron gas has been calculated with 
fixed-node diffusion Monte Carlo, including backflow correlations, for a wide range of electron den- 
sities as a function of spin polarization. We give a simple analytic representation of the correlation 
energy which fits the density and polarization dependence of the simulation data and includes sev- 
eral known high- and low-density limits. This parametrization provides a reliable local spin density 
energy functional for two-dimensional systems and an estimate for the spin susceptibility. Within 
the proposed model for the correlation energy, a weakly first-order polarization transition occurs 
shortly before Wigner crystallization as the density is lowered. 
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The two-dimensional electron gas (2DEG), realized in 
semiconductor heterostructures |l| , exhibits an extremely 
rich phenomenology at low density, where correlations 
play an important role. Being outside the reach of pertur- 
bative approaches, many crucial aspects of this interest- 
ing physics still lack a satisfactory explanation 0] • Under 
these circumstances, valuable information can be gained 
from simplified models, like the ideal 2DEG (strictly 2D 
electrons interacting via a 1/r potential within an uni- 
form, rigid neutralizing background). At zero tempera- 
ture, the state of this system is entirely specified by the 
coupling parameter r s = 1/ y/nna b , where n is the den- 
sity and as the Bohr radius, and the spin polarization 
£. Even for such a simple model the theory has greatly 
benefitted from numerical work, resulting -at least for a 
limited range of physical properties and/or system sizes- 
in benchmark results, input quantities for approximate 
theories, and aspects of the ground-state phase diagram. 
The Quantum Monte Carlo (QMC) method Q, for exam- 
ple, has accurately predicted the critical density for 
Wigner crystallization of the ideal 2DEG,later observed 
experimentally Q , and has characterized the polariza- 
tion transition as being weakly first-order and occurring 
shortly before crystallization as the density is lowered. 

In this work we present extensive QMC simulations of 
the liquid phase within the whole range of density and 
polarization, and provide an analytic expression for the 
correlation energy e c as a function of r s and £. Previous 
estimates of the £ dependence were based on interpola- 
tion conjectures between the energy of the paramagnetic 
and of the fully polarized liquid 0, but, as shown 
later, these estimates may significantly depart from the 
calculated energy at intermediate polarization. 

The interest in the C dependence stems from its use 
in the spin-density functional theory (SDFT) of two- 
dimensional systems 0, 0, in the development of 
density functionals in the presence of magnetic fields ^lj > 
and, more generally, in the study of the spin-polarization 
transition, whose role on transport properties 0] is 
prompting intense experimental investigation^^- 

Our combination of numerical results and known an- 



alytic properties yields an estimate of the paramagnetic 
susceptibility. At low density, this is an extremely diffi- 
cult quantity to evaluate with approximate theories [l^j ] 
due to the need of disentangling the tiny -but important- 
effect of quantum statistics from the huge effect of cor- 
relation. Finall y, s ince we include the effect of back- 
flow correlations [l^ as explained below, our e c (r s , C) pro- 
vides a more accurate phase diagram than previously re- 
ported 0, HI • Hartree atomic units are used throughout 
this work. 

Numerical results - Our calculations use standard 
fixed-node diffusion Monte Carlo (FN-DMC) @], which 
projects the lowest-energy eigenstate $ of the many-body 
Hamiltonian with the boundary condition that $ van- 
ishes at the nodes of a trial function '5. The algorithm 
simulates the imaginary-time evolution by a branching 
random walk of M copies of the iV-electron system, using 
a short-time approximation of the importance-sampled 
Green's function 3] . For each of the densities correspond- 
ing t o Ts = 1) 2, 5, 10 we have considered about 20 values 
of TvHa and 10-12 polarizations ( = (N^ - N l )/N. A 
few simulations have been done also for r s =40, ( = 1. 
The electrons are placed in a square box with periodic 
boundary conditions. Long-ranged interactions are dealt 
with by a model potential (Eq. (68) of Ref. [lfij |) which 
has been shown to give smaller finite-size corrections than 
Ewald sums for the electron gas. The bias introduced by 
the finite population M of walkers and the finite time step 
t was evaluated for selected systems and interpolated for 
all the others, according to standard procedures. To es- 
timate the difference A between the energy eAr(r s ,£) of 
the finite system and its thermodynamic limit e(r s , £), de- 
fined as TV — > oo at fixed density, we adopted a less usual 
strategy. Rather than a separate size extrapolation for 
each density based on variational energies 0, 0, Q] , we 
performed a global fit directly based on FN-DMC ener- 
gies, which exploits two physically motivated ingredients: 
(i) the Fermi-liquid-like size correction 01 

A(r s ,C,iV)= ejv(r SJ C)-e(r SJ C) = 

6(1 + AC 2 )[ t N (r Sl C) - i.(r„ C)] ~ (v + lC 2 )/N{\) 



2 



[tjsr is the energy per particle of TV noninteract 
trons confined in a 2D box of side L = (N/n) 
periodic boundary conditions, t s = (1 + ( 2 )/2r 2 is equal 
to limjv^oo% with n kept fixed, and S, A, 77,7 are r s - 
dependent parameters]; (ii) an analytic expression for 
e(r s ,£), detailed below, which involves 12 more free pa- 
rameters. On these grounds all the FN-DMC energies 
£Ar(r s ,C) calculated here for r s = 1 to 10, plus those 
of Ref. |3| for r s = 20 and 30, plus the energy (inde- 
pendently extrapolated to the thermodynamic limit) for 
r s = 40, C = 1 -a, total of 122 data- formed the input 
for a best fit of the 36 free parameters, 24 of which dis- 
appear from the final analytic expression since they only 
concern the size extrapolation. This fit yields a reduced 
X 2 of 3.8. More details will be reported elsewhere. 

The fixed-node approximation is variational in charac- 
ter |3j, and its accuracy depends on the nodal structure 
of 'J. We choose a Slater- Jastrow trial function ^f(R) = 
J(R)D^(R)D L (R) where ft={ri, . . . ,r N } represents the 
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electronic coordinates, J(R) — exp Y^i j u si,s :i (i"ij) 
symmetric Jastrow factor with differe nt p seudopotentials 
u Si , S j for like- and unlike-spin pairs |l8|, and Z)T(i) is a 
Slater determinant for spin-up(down) electrons. 

The standard choice with homogeneous systems is to 
use plane waves (PW) as single-particle orbitals: -Dp^ = 
det [expfzk, ■ r-,-1][l5j. However, within the fixed-node 
approximation, better results are obtained with back- 



flow (BF) correlations in the wave function [l4j , D 



tU) 

BF — 

det [exp(zk. J • x,-)], with x, : = r 4 + Vbf (nj ) fa - r,-)- 
Since simulations with the BF wave function are con- 
siderably more demanding than with PW determinants, 
we calculated BF energies only for ( — Q, N — 58 and 
£=1,TV = 57 for each density, including correlation func- 
tions tibf(t) as described in Ref. Jl4j optimized with vari- 
ational Monte Carlo simulations [3j in each case. The ex- 
pected error with BF nodes is much smaller than the dif- 
ference between the PW and BF energies. For r s — 1, an 
exact calculation^!^] shows agreement with the BF result 
within the statistical error. The BF results are compared 
with PW energies in Table |U For other values of TV and 
£ the effect of backflow is estimated as a quadratic inter- 
polation in £ and appended to PW energies, under the 
further assumption that the size dependence be the same 
for BF and PW: the - 120 FN-DMC data used in the 
fit of Eq. (JJJ are the calculated PW energies, corrected 
for the time-step and population bias, and shifted by the 
estimated backflow effect -the infinite-size extrapolation 
being given as a result of the fit. 

As expected |7J , BF correlations lower the energy more 
in the paramagnetic than in the polarized phase. At 
large r s this relative gain of the paramagnetic phase is 
a significant fraction of the difference in energy between 
the two phases. 

Analytic representation for e c (r s , £) - We parametrize 
the energy e(r s ,Q as follows. We first noticed that the 
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-0.2013(1) 
-0.255802(4) 
-0.149134(9) 
-0.0852706(4) 
-0.046241(1) 
-0.031923(1) 



0.13147(2) 
-0.193349(1) 
-0.143520(5) 
-0.084555(2) 
-0.0462385(6) 
-0.0319298(6) 



-0.20372(4) 
-0.25721(3) 
-0.149518(9) 
-0.085427(6) 
-0.046283(1) 
-0.031941(2) 



0.13109(4) 
-0.19359(2) 
-0.143610(7) 
-0.084584(2) 
-0.0462488(8) 
-0.031938(1) 



TABLE I: FN-DMC energies with plane wave or backflow 
nodes. Data pertain to simulations with TV = 58 for £ = 
and N = 57 for £ = 1, M = 200 and r = 0.002, 0.01, 0.1, 0.3, 
1.0, 2.0 in order of increasing r s . 
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-0.1925* 
0.0863136* 

0.0572384 
1.0022 

-0.02069 

0.33997 
1.747 x 10 -5 



0.117331* 0.0234188* 

-3.394 x 10" 2 -0.037093* 

-7.66765 x 10" 3 * 0.0163618* 
0.4133 1.424301 
0* 0* 
6.68467 x 10" 2 0* 



7.799 x 10" 



1.163099 



TABLE II: Optimal fit parameters for the correlation energy, 
as parametrized in Eqs. @ and @. Values labelled with * are 
obtained from exact conditions. The parameter D; = —AiHi 
is not listed (see text); the parameter j3 is equal to 1.3386. 



spin-polarization dependence of the exchange-correlation 
energy e xc = e—t s , as given by the DMC data, is very well 
described by the simple form cq + ci( 2 + C2C 4 for r. > 5. 
On the other hand, the known high-density limit (2fJ,|2l|, 



e X c(r s ,Q = e x {r s ,Q + a (Q + b (C)r s lnr s + 0(r s ), (2) 

contains non-negligible contributions from higher pow- 
ers of (: the dominating exchange term e x is equal to 
-2\/2[(l + C) 3/2 + (1 - C) 3/2 ]/37rrv Since we want to 
interpolate the energy between high and low density, 
we choose a functional form for the correlation energy 
(c = (xc — £x which quenches the contributions from e x 
beyond fourth order in ( as r s increases, 

e c (r s C) = (e-^«-l)4 6) fa,C) 

+ a {r s ) + ai (r s )( 2 +a 2 (r s )(\ (3) 

where 4 6) (r s ,C) = e«(r„C) - (1 + |C 2 + ^M^.Q) 
is the Taylor expansion of e x beyond fourth order in (. 
Since the first term in the rhs of Eq. © contains pow- 
ers 6 and higher of Cj it is immediate to identify the 
function ao(r s ) as the correlation energy at zero po- 
larization, ao(r s ) = e c (r s ,0). Furthermore, 2ai(r s ) = 
[<9 2 e c (r s ,C)/<9C 2 ] c=0 (spin stiffness), and 24a 2 (r s ) = 
[<9 4 e c (r s , C)/9C 4 ]^_ • Our choice for the functions aj(r s ) 
is a generalization of the Perdew-Wang [22j form to the 
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FIG. 1: Left panel: (^-dependence of the correlation energy 
in the high-density limit. The result of the present work 
is compared with the exact values from Ref. |2fl | and with 
the exchange-like interpolation of Refs. ||. Right panel: £- 
dependence of the total energy at the polarization transition 
density, r s = 25.56. The value at £ = has been subtracted, 
and the result multiplied by 10 5 . Our result is compared with 
the exchange-like interpolation. 
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which features both the subleading contributions of the 
expansion (J2J and terms in rj 1 and r s 3 ^ 2 for r s — ► oo|23|. 
With suitable constraints, which leave only 12 indepen- 
dent parameters in Eqs. (J3J and J2J, our correlation en- 
ergy satisfies exactly several known high-density and low- 
density limits: (i) the exact values |2C1 121| of ao(C) and 
&o(C) at C — and £ = 1 in the small-r s expansion of 
e c shown in Eq. J3J); (ii) the vanishing of the correlation 
energy in the r s — > oo limit, which simply implies that 
Di = —AiHi] (in) the requirement that e(r s ,() be inde- 
pendent of £ for r s — > oo. The latter condition is enforced 
up to order 0(rj 2 ), thus recovering the low-density be- 
havior e — * —m/r s + njr^J 2 + 0(rj 2 )^^ with positive 
m and n independent of £. 

The optimal values of the parameters are listed in Ta- 
ble [n] In the high-density limit, we can compare the 
^-dependence of our correlation energy with the exact 
one |2fJ|, finding very good agreement, as shown in the 
left panel of Fig. [IJ where the widely-used exchange- 
like interpolation [8j is also reported. In the right panel 
we instead see the ^-dependence of the total energy at 
r s = 25.56 where, according to our results, the transition 
to the fully polarized gas occurs. At such low densities, 
the exchange-like interpolation scheme (here performed 
using our energy values at £ = and C — 1) significantly 
deviates from the QMC result. Both predict a sudden 
transition from the ( = to the £ = 1 fluid, but the 
energy barrier between the two phases given by QMC is 
more then an order of magnitude smaller, which reflects 
in a very large value of the spin susceptibility. 




FIG. 2: Spin-susceptibility for the 2DEG as a function of r 3 . 
The present result is compared with the exchange-like approx- 
imation 8] , with two quadratic interpolations . J| based on the 
total energy (QI etot ) and on the correlation energy (QI ec ), and 
with the Yarlagadda and Giuliani (YG) calculation. 



Spin susceptibility - In our parametrization, the spin 
susceptibility \ of the ideal 2DEG is simply: 
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where a.\ is given by Eq. (0} and Table ITU As mentioned, 
it turns out that \ is a very delicate quantity, so that 
different theoretical predictions significantly depart from 
each other at r s as low as 3 or 4 [l3j. 

In Fig. |3 we compare our result with other estimates 
of x- The result of the Yarlagadda and Giuliani (YG) 
calculation blows up already at r s ~ 4. The ten- 
dency to predict a polarization transition at too high 
densities is shared by several approximate theories |l3|. 
The exchange-like interpolation [8| is significantly lower 
than the QMC estimate, as expected from Fig. 2] Out 
of two recipes given by Tanatar and Ceperley (TC) 0], 
based on a quadratic interpolation of the C-dependence 
of the total energy (QE tot ) and of the correlation energy 
(QI £c ) respectively, the former (which yields by defini- 
tion a divergent x at the polarization density) is quite 
accurate for r s < 10, whereas the latter greatly under- 
estimates the spin-susceptibility at all but the smallest 
r s (here the quadratic interpolation has been performed 
using our energy values at ( = and ( = 1). 

Our spin susceptibility is related to the second deriva- 
tive of the model correlation energy, which incorporates 
an optimal interpolation of the QMC results and is con- 
strained by known limiting behaviors at very low and 
very high densities. For this reason it represents a consid- 
erable progress over existing theories, providing a sound 
reference for further studies. 

Nevertheless, the precise value of x at very large r s 
(say, above 20) has to be taken with some caution. When 
the C-dependence of e(r s ,() is extremely weak (see the 
lower curve in the right panel of Fig. ^) the calculation 
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FIG. 3: r s -dependence of the ground-state energy e of the 
2DEG for the paramagnetic (£ = 0) and ferromagnetic (£ = 1) 
fluid phases, and for the Wigner crystal as given in Ref. 
€m = — 1.1061/r s is the Madelung energy. 



of its derivatives is beyond the accuracy of the present 
calculation, due to the fixed-node bias, the assumption of 
quadratic contributions from BF correlation, possible un- 
certainties from the size extrapolation, statistical noise, 
and the chosen functional form of e c , Eq. Jjy|. The very 
nature of the polarization transition, wealky first-order 
according to Fig. is clearly based on the above approx- 
imations and assumptions. 

Phase diagram - The above results allow us to draw the 
zero-temperature phase diagram shown in Fig. [21 rela- 
tive to the Wigner crystal (WC), the paramagnetic liquid 
(PL) and the ferromagnetic liquid (FL). Intermediate- 
polarization phases never represent the stable phase in 
2D (Ref. [7j and Fig. . For the fluid phases, the energy 
e — t s + e xc is given by Eq. Q with the parameters of 
Table|nJ The energy of the crystal is taken, instead, from 
Ref. [5j , since neither backflow nor spin polarization play 
a significant role in the solid phase |2J] . 

Similar FN-DMC studies, using plane- wave nodes, 
have been previously performed |5j. According to 
Ref. [4j, crystallization occurs directly from the PL, 
although, at freezing, the energies of all three phases 
are very close to each other. Subsequent PW-based 
simulations^ Q revised slightly upwards the energy of 
the PL, but this was enough to alter the previous result, 
and to predict a small stability window for the FL 5]. 
Our backflow calculations lower back the energy of the 
PL relative to the FL, with the effect of shrinking, but 
apparently not eliminating, the density range where the 
FL phase is stable: the ideal 2DEG undergoes a polar- 
ization transition at r s ~ 26, and the polarized liquid 
crystallizes at r s ~ 35. 
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